function [VT, VE, ring1, areas, bnd_flag] = build_geom_info(V,T,E,ET)
% function [VT VE, ring1, areas, bnd_flag] = build_geom_info(V,T,E,ET)
%

nt = size(T,1);
nv = size(V,1);

%% 1. build VT first
i = zeros(3*nt,1);
j = i;
i(1:nt) = T(:,1);
j(1:nt) = 1:nt;
i(nt+1:2*nt) = T(:,2);
j(nt+1:2*nt) = 1:nt;
i(2*nt+1:3*nt) = T(:,3);
j(2*nt+1:3*nt) = 1:nt;
VT = sparse(i,j,ones(3*nt,1),nv,nt);

%% 2. find the boundary Nodes
bdr = find(ET(:,2)==0);
bnd_flag = zeros(nv,1);
bnd_flag(E(bdr,1)) = 1;
bnd_flag(E(bdr,2)) = 1;


%% 3.build first ring
ne = size(E,1);
i = zeros(2*ne,1); j = i;
i(1:ne) = E(:,1); i(ne+1:2*ne) = E(:,2);
j(1:ne) = E(:,2); j(ne+1:2*ne) = E(:,1);
ring1 = sparse(i,j,ones(2*ne,1),nv,nv);

%% 4. calculate areas for each triangle
x21 = V(T(:,2),1)-V(T(:,1),1);
y21 = V(T(:,2),2)-V(T(:,1),2);
x31 = V(T(:,3),1)-V(T(:,1),1);
y31 = V(T(:,3),2)-V(T(:,1),2);
areas = 0.5*(x21.*y31 - x31.*y21);

%% 5.% at last build matrix VE, which indicate the edges around some vertex
%% 
 VE = comput_VE(size(V,1),E);
 
end


%%%%
function VE = comput_VE(nV,E)
    nE = size(E,1); 
    i(1:2:2*nE-1) = 1:nE;
    i(2:2:2*nE)   = 1:nE;
    j(1:2:2*nE-1) = E(:,1);
    j(2:2:2*nE)   = E(:,2);
    s = ones(2*nE,1);
    EV = sparse(i,j,s,nE,nV);

    count = sum(EV)';
    nmaxedge = max(count);

    % begin to build VE
    VE = zeros(nV,nmaxedge);
    for k = 1:nE
        v1 = E(k,1);
        idx = VE(v1,1);
        VE(v1,1) = idx + 1;
        VE(v1,idx+2) = k;

        v2 = E(k,2);
        idx = VE(v2,1);
        VE(v2,1) = idx + 1;
        VE(v2,idx+2) = k;
    end
    if (sum(count - VE(:,1)) ~= 0)
        fprintf('error in mesh data VE!\n');
    end

    % sort VE in an anti-clock term
    % for k = 1:nV
    %     nedge = EV(k,1);
    %     edges = EV(k,2:nedge+1);    
    %     idx = sort_edge(V,T,E,ET,edges,k);
    %     EV(k,2:nedge+1) = edges(idx);  % reorder the edges' index
    % end

end